home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / wathen_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  2.1 KB  |  86 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Wathen matrix - a finite element matrix (random entries).
  4.  
  5. // Syntax:    A = wathen ( NX , NY )
  6.  
  7. // Description:
  8.  
  9. //    A is a sparse random N-by-N finite element matrix where 
  10. //    N = 3*NX*NY + 2*NX + 2*NY + 1. 
  11.  
  12. //    A is precisely the `consistent mass matrix' for a regular
  13. //    NX-by-NY grid of 8-node (serendipity) elements in 2 space
  14. //    dimensions. A is symmetric positive definite for any
  15. //    (positive) values of the `density', RHO(NX,NY), which is
  16. //    chosen randomly in this routine. In particular, if 
  17. //    D = DIAG(DIAG(A)), then  
  18.  
  19. //              0.25 <= EIG(INV(D)*A) <= 4.5
  20.  
  21. //    for any positive integers NX and NY and any densities
  22. //    RHO(NX,NY). This diagonally scaled matrix is returned by
  23. //    WATHEN(NX,NY,1). 
  24.  
  25. //      Reference:
  26. //        A.J. Wathen, Realistic eigenvalue bounds for the Galerkin
  27. //        mass matrix, IMA J. Numer. Anal., 7 (1987), pp. 449-457.
  28.  
  29. //    This file is a translation of wathen.m from version 2.0 of
  30. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  31. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  32.  
  33. //-------------------------------------------------------------------//
  34.  
  35. wathen = function ( nx , ny , k )
  36. {
  37.   local ( nx , ny , k )
  38.  
  39.   if (!exist (ny)) 
  40.   { 
  41.     error("Two dimensioning arguments must be specified.");
  42.   }
  43.   if (!exist (k)) { k = 0; }
  44.  
  45.   e1 = [6, -6, 2, -8; -6, 32, -6, 20; 2, -6, 6, -6; -8, 20, -6, 32];
  46.   e2 = [3, -8, 2, -6; -8, 16, -8, 20; 2, -8, 3, -8; -6, 20, -8, 16];
  47.   e = [e1, e2; e2', e1]/45;
  48.   n = 3*nx*ny+2*nx+2*ny+1;
  49.   A = zeros(n,n);
  50.  
  51.   RHO = 100*rand(nx,ny);
  52.  
  53.   for (j in 1:ny)
  54.   {
  55.     for (i in 1:nx)
  56.     {
  57.       nn[1] = 3*j*nx+2*i+2*j+1;
  58.       nn[2] = nn[1]-1;
  59.       nn[3] = nn[2]-1;
  60.       nn[4] = (3*j-1)*nx+2*j+i-1;
  61.       nn[5] = 3*(j-1)*nx+2*i+2*j-3;
  62.       nn[6] = nn[5]+1;
  63.       nn[7] = nn[6]+1;
  64.       nn[8] = nn[4]+1;
  65.  
  66.       em = e*RHO[i;j];
  67.  
  68.       for (krow in 1:8)
  69.       {
  70.         for (kcol in 1:8)
  71.         {
  72.           A[nn[krow];nn[kcol]] = A[nn[krow];nn[kcol]]+em[krow;kcol];
  73.         }
  74.       }
  75.  
  76.     }
  77.   }
  78.  
  79.   if (k == 1)
  80.   {
  81.     A = diag(diag(A)) \ A;
  82.   }
  83.  
  84.   return A;
  85. };
  86.